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Abstract 

The Hierarchical Reference Theory (hrt) and the Self-Consistent Ornstein- 
Zernike Approximation (scoza) are two liquid state theories that both fur- 
nish a largely satisfactory description of the critical region as well as phase 
coexistence and the equation of state in general. Furthermore, there are a 
number of similarities that suggest the possibility of a unification of both 
theories. As a first step towards this goal we consider the problem of com- 
bining the lowest order 7 expansion result for the incorporation of a Fourier 
component of the interaction with the requirement of consistency between 
internal and free energies, leaving aside the compressibility relation. For 
simplicity we restrict ourselves to a simplified lattice gas that is expected 
to display the same qualitative behavior as more elaborate models. It turns 
out that the analytically tractable Mean Spherical Approximation is a so- 
lution to this problem, as are several of its generalizations. Analysis of the 
characteristic equations shows the potential for a practical scheme and yields 
necessary conditions any closure to the Ornstein Zernike relation must fulfill 
for the consistency problem to be well posed and to have a unique differen- 
tiable solution. These criteria are expected to remain valid for more general 
discrete and continuous systems, even if consistency with the compressibility 
route is also enforced where possible explicit solutions will require numerical 
evaluations. 

1 Introduction 

Both the Self-Consistent Ornstein-Zernike Approximation (scoza, [1, 2, 3]) and 
the Hierarchical Reference Theory (hrt, [4]) have been found to give very accurate 
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results for fluids in thermal equilibrium. In particular, the respective non-linear 
parabolic partial differential equations (PDEs) remain practical in the critical re- 
gion, and their solution gives non-classical, and partly Ising-like, critical indices. 
The PDEs themselves are derived by obtaining the equation of state in two inde- 
pendent ways and using thermodynamic consistency to fix a free parameter in the 
direct correlation function. 

Although both approaches thus appear similar in a number of aspects, there 
are also marked differences: While both make use of the compressibility route 
to thermodynamics, SCOZA combines it with the energy route expression for the 
internal energy, whereas hrt, inspired by momentum-space renormalization group 
theory, relies on what might be called the fluctuation route to the free energy 
instead. In addition to the density p, the independent variables in the PDEs are 
therefore the inverse temperature (3 = 1/ksT and the momentum-space cutoff Q, 
respectively. Starting from a reference system of known properties at vanishing (3 
or at high Q, the attractive interaction is then turned on by gradually increasing 
its strength (3 (scoza), or by including its Fourier component of wavenumber Q 
at constant temperature until the full system is recovered in the limit Q — > 0. 

In the present contribution we want to investigate the possibility of combin- 
ing both approaches by imposing thermodynamic consistency of the internal and 
free energies to obtain a differential formulation where both the strength of the 
interaction and its Fourier components are added successively. This introduces an 
additional constraint so that a further free parameter in the closure relation can be 
determined. In combination with the compressibility route there are then two free 
parameters that may be used, e. g., to describe both the range and the amplitude 
of a contribution to the direct correlation function c(r); it might then be possi- 
ble to describe the long-ranged tail that c(r) is known to develop in the critical 
region. However, our preliminary investigations point to several difficulties with 
this approach near phase coexistence. The crucial problem is the vastly different 
behavior of hrt and SCOZA in this part of the phase diagram: hrt has a solution 
at all densities, but spinodal and binodal coincide so that the two phases coexist 
at infinite compressibility; SCOZA, on the other hand, gives distinct binodal and 
spinodal curves but does not have a solution inside the spinodal. At this point it 
is unclear how these differences might be resolved by a simple modification of the 
direct correlation function. Even at temperatures above the critical one, (3 < f3 c , 
where full self-consistency should be attainable, we will be confronted with dif- 
ferential equations that are of first order in temperature and cutoff and of second 
order in the density. Numerical solution of such a pde can be demanding. 

To simplify our task and gain some insight we here consider a simplified lattice 
gas model, and we limit ourselves to combining the energy and fluctuation routes 
without using the compressibility at all. Then there is only one free parameter left 
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that can be determined as a function of (3 and Q. After introducing our model and 
establishing the pde (section 2), we show that consistency is achieved with the 
mean spherical approximation (msa) and several generalizations of it (section 3). 
The solution can therefore be given in closed analytic form, but the price to pay for 
omission of the compressibility is the presence of the MSA features, viz., "irregular" 
isotherms in the critical region and lack of a well-defined critical behavior for the 
lattice gas case [5]. 

In general, an analytical solution is also unlikely to be known beforehand if 
other closures are used. In that case one should seek the solution numerically by 
integrating the equations for the characteristic curves of the pde. Especially when 
the compressibility relation not considered here is incorporated, a numerical solu- 
tion will be required, and we expect the properties of the characteristics analyzed 
here to be useful more generally. After establishing relations and implications for 
more general situations (section 4.1) we once more turn to the MSA equations. 
It is found that the pde admits a more general solution. Usually, the reference 
system also does not determine the solution at Q = for finite temperature due to 
the behavior of the characteristics (section 4.2). A generalized version of the pde 
that does not make any assumptions on the form of the direct correlation function 
of the system at cutoff Q and inverse temperature (5 other than that it depends 
on some unknown parameter function A(/3, Q) is considered in section 4.3. It is 
then found that these deficiencies can be peculiar to the specific closure employed. 
Based on these findings, we finally state and consider two complementary neces- 
sary conditions a parameterization of the solution must fulfill for the consistency 
problem to be well posed and physically reasonable (section 5). These conditions 
have the advantage that they can be checked in the high temperature limit where 
significant simplifications are possible. They are also sufficiently general to apply 
to more realistic model systems, both discrete and continuous. In view of the re- 
lation between the properties of a non-linear diffusion equation with two time-like 
parameters on the one hand and the first-order pde we consider in the present con- 
tribution on the other hand, the criteria found for the selection of suitable closure 
relations remain relevant even when the compressibility route to thermodynamics 
is also taken into account. 

2 Adaptation of theory and model 
2.1 Model 

In the interest of simplicity, we here consider only a lattice gas (or Ising model). 
The main advantage of this choice is that the core condition of vanishing pair dis- 
tribution function g(r) = h(r) + 1 affects only the single point at r = ]f| = 0. In 
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addition, we assume a sufficiently long-ranged interaction so that the anisotropy 
imposed by the lattice structure can be neglected. Besides, an even more long- 
ranged tail appears in the interaction for intermediate Q, cf. [6] and appendix D.l 
of [7]. When following the hrt recipe of successively including Fourier compo- 
nents of the interaction, we will therefore do so without regard to the geometry 
of the Brillouin zone. In the same spirit we will often speak of finite or infinite 
cutoffs and wavenumbers instead of specifying the corresponding surfaces within 
the Brillouin zone explicitly [8, 9], and we will also write Qoo for the maximum 
cutoff in the calculation. It is not expected that these simplifications should af- 
fect the results qualitatively. Furthermore, the main conclusions of section 5 are 
manifestly independent of these assumptions. 

The system is assumed to differ by a perturbing attractive interaction —ip(r), 
ip(r) > 0, from a reference system of known properties. In the lattice case the 
simple hard core lattice gas can serve as the reference system. We can then use 
the inverse range 7 of the attractive interaction as a perturbing parameter [10, 11]. 
The zeroth order contribution past the reference system is the mean-field attractive 
van der Waals term, corresponding to a pressure contribution of — \ tp(0) p 2 , where 
a tilde indicates Fourier transformed quantities. In hrt this corresponds to the 
zero-loop diagram [4]. To next order in 7 Hemmer found [10] 



/ = -C J In (1- fl(k)v(k)) d 3 A;oc7 3 (1) 



with 



2 \2tc 
P = ^T[ 

Here, I is —f3 times the first order contribution to the Helmholtz free energy per 
unit volume, ks is Boltzmann's constant, and T is the temperature. The integral 
is to be extended over the Brillouin zone of the lattice that we here assume to be 
the simple cubic one for which C J d 3 k = \. 

The function p{r) is related to the direct and total correlation functions, co(r) 
and ho(r), of the reference system (indicated by the subscript 0) by 

fl(k) = p + p 2 h (k) 

= P (2) 
1 - pc (k) 

Here we have used the Ornstein-Zernike equation [12] 

h(r) = c(r) + p J h{f — f ') c(f') d 3 r', 
(l + ph(k))(l-pc(k)) = l. 



(3) 
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In the lattice case, the integral is to be replaced with a sum. 

In the limit 7 — » 0, the function v(k) usually reduces to a narrow peak of width 
7 at k = 0. In this limit the expression (1) is the exact correction to the mean 
field result. In hrt this corresponds to the infinite sum of one-loop diagrams [4]. 

The interaction v(k) may also have a narrow peak of width 7 at some other 
position, say, around k = Q. In r space this corresponds to an infinitely weak 
but oscillating interaction with wave vector Q, and eq. (1) is still exact in the 
limit 7^0. On the other hand, the renormalization procedure of hrt consists 
of adding narrow pieces of ip of width — dQ > around shorter and shorter wave- 
vectors Q at constant temperature. With d 3 k = An k 2 dk and after division by dQ, 
eq. (1) then becomes 

81 

— = AnCQ 2 - j,(Q)v(Q)) ■ (4) 

This is also the first equation of the hrt hierarchy, provided \x is derived from the 
structure of the system at slightly higher cutoff Q + \dQ\ rather than from the 
reference system at Q^. Eq. (4) is then formally exact, but it does not specify 
how p, changes as Fourier components of ip are added and the cutoff Q approaches 
zero. The formal answer to this problem is provided by the higher equations of 
the HRT hierarchy or some other expansion for \x. In practical applications of 
hrt only eq. (4) is used, and the evolution of \x is determined by introducing a 
free parameter into the correlation functions and fixing it by the requirement of 
consistency with the compressibility route. 

One can also make a small change d(3 > in the inverse temperature and so 
effectively increase the interaction by an amount —dv(k) = ip(k) df3. Taking the 
limit dp — > we obtain 

pUl -dp~ c J i-mm { ) 

The quantity u\ is nothing but the configurational internal energy per particle 
beyond the zeroth order mean field term (besides a self-energy term — 1^(0) in- 
cluded here from the p-term in /i) due to a structure or correlation functions given 
by 

With n given by eq. (2) this is exact only in the limit P — > where the reference 
system is recovered. At higher P, it is only the first order correction in 7 for the 
long-range part of the pair correlation function [10, 11]. To higher order, or at 
finite P, it is again not obvious how p,(k) changes with temperature. Here the 
SCOZA recipe is to introduce some parameter into the closure that is also fixed by 
the requirement of consistency with the compressibility route. 
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2.2 Specialization to lattice gas 

In the case of a lattice gas or the Ising model, the above expressions are simpler. 
The direct correlation function for the reference system vanishes except at r = 
so that jl(k) is constant, jl(k) = /i; and from the pressure p = — ln(l — p)/ (3 of the 
hard core lattice gas serving as reference system we obtain its value as 

^ = #°) = Q fr) 1 = p^-p)- (7) 

As a simple approximation we can now replace the \i of eq. (7) by some effective 
value fi e that may be determined by imposing thermodynamic consistency between 
eqs. (4) and (5); in general fi e will be a function of (3 and Q. Such a change may 
cause the core condition h(0) = —1 to be violated, which might be corrected by 
the introduction of an additional parameter. For the qualitative features we are 
interested in here, however, the core condition is not expected to be important, 
and we do not impose it below [13]. 

As Fourier components are added to the interaction at successively smaller 
cutoff Q, the interaction vanishes for k < Q (except for the mean field term at 
k — 0), and the integral in eq. (5) must be restricted to that part of the Brillouin 
zone where k > Q. Furthermore, a small change dQ in the cutoff corresponds to the 
addition of a weak long-ranged oscillating tail in r space, which also contributes 
at r = 0. By setting the potential ip at r = equal to zero the corresponding 
unphysical contributions to the internal energy from within the hard core can be 
avoided approximately, and exactly in the limits Q — > Qoo and Q — > 0, even when 
the core condition is not fulfilled. With these considerations and using spherical 
symmetry eqs. (4) and (5) become 

— =inCQ 2 ln(l-z4:(Q)) 

m cr z _m < 8 > 

d(3 PJk>Ql-z^(k) 

with 

z = pe!3. (9) 

The above equations define a pde for p e ((3, Q) that is to be solved for all Q > 0, 
/3 > 0. In particular, the solution should give the structural and thermodynamic 
properties of the target system at Q = and all (3 > 0. In order to turn this 
into a well-defined problem we also have to impose some boundary conditions. 
The natural choice is to demand that the reference system (/x e = /i) should be 
recovered at f3 — and at Q — Q^. Whether this is sufficient to determine the 
solution throughout the domain of the PDE, however, is not obvious and depends 
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not only on the general consistency problem but also on the particular closure (the 
replacement of \x by p e ) and the toy model chosen. We will return to this point in 
greater detail in section 4. 

At first sight, an analytic solution of the given equations is not obvious. But as 
we will show in the following section, the MSA solution and generalizations thereof 
indeed fulfill eq. (8). On the other hand, the multiplicity of solutions found already 
indicates that the reference system boundary conditions (that are fulfilled in all 
cases) do not determine a unique solution, cf. sections 3.2 and 4.2. 



3 MSA and MSA-like solutions 



3.1 Specific solution 

The internal energy corresponding to the MSA correlation function can be inte- 
grated explicitly to yield the free energy by utilizing general MSA expressions [3]. 
The result is 

I = -4nc[ ln(l-zi>(k)) k 2 dk + J -\ ln(l + 2J). (10) 
Jk>Q v ' 2 

Here J is the integral appearing in eq. (8), 



_zifj{k) ( U ) 

and the core condition determines z as 



4vrC f Z ^Z> k 2 dk, 
Jk>Q 1 - zip(k) 



z = 



1 + 2J (12) 
fi =p(l- p ). 

By differentiation of eq. (10) with respect to (5 and Q and using eqs. (11) and (12) 
one finds that eq. (8) is fulfilled. 

The pair correlation function consistent with solution (10)— (12) is obtained 
from eq. (6) by replacing fl{k) by /i e = z/j3 and restricting the interaction to 

k>Q, 

p + p 2 h(k) - 

i-zi/>(k)e(k-Q) 

where is Heaviside's function. Integrating the above over the Brillouin zone and 
using eq. (11) one easily obtains 

p + p 2 h(0) = ^(l + 2J)=p = p-p 2 , 
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which shows that the core condition is fulfilled. 



K = -=Cl : ^L d^ 



3.2 More general solution 

The MSA solution just discussed is one specific solution of eq. (8). A more general 
solution is suggested by previous work on the resummation of 7 ordering for fluids 
[5, 14], a procedure that transforms the serious divergence at the mean field critical 
point into MSA type criticality while maintaining self-consistency between the free 
and internal energies. 

The generalized solution is then 

/oo 
In (1 - fi e v(k)) d 3 k - A n K n+1 (13) 

with (z = f3fi e ) 

= £ = c / T 

He Jk>Q 1 

00 

fi e = H + A n Rn 
n=l 

Here the coefficients A n are independent of f3 and Q and implicitly define /z e . The 
A n must be fixed by some suitable boundary conditions. These must be imposed at 
finite temperature and cutoff as K vanishes at f3 — and at Q — Qoo- This shows 
that the reference system alone does not define a unique solution in the present 
case, cf. end of section 2.2. The MSA solution (10) corresponds to the more general 
eq. (13) with one specific choice of the A n . By differentiation it is easily verified 
that the above expression for /, like expression (10), solves both equations (8). 

3.3 Wavenumber dependent direct correlation function 

Solution (13) can easily be extended to some very specific k dependent /i e and z. 
Replacing /i e by jle(k), the Fourier transform of /i e (r) that now extends to outside 
the core, we can write 

00 

f , e (r)= f ,5(r) + J2Mr)K(r) n . (15) 

71=1 

This corresponds to graphs with n parallel K(r) bonds. K(r) is the chain bond 
with renormalized hypervertex /x e (r), and its Fourier transform is 

k ^ = \ i-l%m * {k - Q) ' (16) 
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the generalization of the integrand in eq. (14). Eq. (13) now changes into 



» oo „ 

I = -C In (1 - jl e (k) v(k)) d 3 A: - V \ A n (r) K(r) n+1 d 3 r, (17) 

Jk>Q n=l U+l J 

and the pde (8) becomes 

81 

— =AnCQ 2 Hl-UQ)v(Q)), 

| = 4,§ / r ^« T ^. (18) 
•9/3 P Jk>Q I - He(k) v(k) 

As mentioned at the beginning of section 2, we neglect the lattice structure for our 
simplified model. Of course, for a true lattice system S(r) in eq. (15) should be a 
Kronecker 5 at the origin, and the spatial integral in eq. (17) should be a lattice 
sum. 

The above solution is suggested by the work of H0ye and Olaussen on the two- 
dimensional Coulomb gas where the well-known Kosterlitz-Thouless phase transi- 
tion was evaluated on the basis of a graph expansion [15, 16]. 

Again it can be verified by differentiation that expression (17) solves eq. (18) 
in the same way expression (13) solves eq. (8). To do so one also needs the identity 



J f(r)g(r)d 3 r = 2C J f(k)g(k) d 3 k. 



4 Characteristics of the differential equations 
4.1 First-order PDEs and characteristics 

So far we have only considered specific cases where exact solutions are known. In 
more general situations the characteristics can be used. These are curves along 
which the pde of interest is equivalent to a set of coupled ordinary differential 
equations. In the case of a first order pde considered here, the solution can always 
be obtained by integrating the characteristic equations, starting from some point 
where the solution is known from some boundary condition. This is exact and 
does not incur any loss of information or generality, and it generates a solution 
that is necessarily differentiable along the characteristic. 1 

Clearly, the solution along a full characteristic depends on only a single point 
at the boundary. This reflects the well-known capacity of first-order PDEs for 
discontinuous solutions. Another way in which discontinuities may arise is from 

1 See textbooks on PDEs such as, e. g., Leon Lapidus, George F. Pinder, Numerical solution 
of partial differential equations in science and engineering, New York (Wiley) 1982. 
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crossings of characteristic curves: the solution then becomes multivalued, and 
a resolution restoring uniqueness naturally leads to a shock where the solution 
is no longer differentiable. A complementary difficulty one may encounter are 
rarefactions, i. e., regions not entered by any characteristics. In this case the 
solution is not defined there, although there is often a natural weak formulation 
of the pde that can be used instead. 

In the present work we consider consistency of the free and internal energies 
at inverse temperature f3 when gradually turning on Fourier components of the 
interaction at ever smaller wavelength Q. Neither shocks nor rarefactions are 
expected on physical grounds, i. e., the solution of the pde should both exist and 
be differentiable for j3 > and Q > 0, at least away from the critical point and 
outside the spinodal. Using subscripts to denote partial derivatives, a pde like 
eq. (8) is naturally written in the form 

V x (x,y) = X(x,y; X(x,y)) 
V y (x,y) = Y(x,y;X(x,y)). 

Here A is an unknown function of the independent variables x and y. As we will 
see, these equations imply a first order pde for A the solution of which also gives 
the quantity of interest, 

The above relations are not in the standard form F(x,y, ty x , fy y ) = for a 
first order PDE. The usual expressions for the characteristics are thus not directly 
applicable. One option is to invert, say, X(x, y\ A) with respect to A. The resulting 
expression for A can then be inserted into the equation for ty y , allowing the usual 
equations for the characteristics of non-linear PDEs to be used. The disadvantage 
of this approach is that it leads to highly involved expressions. 

A simpler way of obtaining the characteristic equations is by cross differentia- 
tion of the PDE: Setting ^ xy = ^> yx gives 

X y + X X X y — Y x + Y x X x . 

This is a quasi-linear first-order PDE for A of the form a X x + b X y = c, for which 
the characteristic equations are dx/a = dy/b = dX/c. In combination with the 
pde for ^ itself we immediately find the set of characteristic equations, 

dty y _ dfr 

Y x ~ X x ~ X y - Y x ~ X x Yx - Xx Y x ~ X y Y x -X x Y y ~ ^ x Y x -^ y X x ' 

(19) 

At every point (x, y) along a characteristic the ratios of these differentials deter- 
mine the direction of its tangent in the (x, y) plane as well as the corresponding 
changes in A, ^ x , ^ y , and * along the curve. The equations can easily be solved 
numerically, e. g., by predictor-corrector methods. In the present contribution, 



dx dy dX d^ 
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however, we will not concern ourselves with numerical evaluations. Suffice it to 
say that a straightforward implementation of these equations may require very 
small step sizes, and that an accurate evaluation of the direction of the tangent 
may be difficult, in the vicinity of points where X x = Y x = 0; discretizations treat- 
ing X x and Y x as perturbations relative to, e. g., X xx and Y xx can be employed 
there. 

An important consequence of the equations given above concerns the orienta- 
tion of the characteristics in the (x, y) plane at special points: Evidently, wherever 
X\ = and Y x ^ 0, the tangent is parallel to the x axis, i. e., y — const. By 
locating the zeros of the A derivatives of the right hand sides of the pde we can 
thus rapidly gain a qualitative overview of the geometry of the field of character- 
istic curves. Considering the path of integration in the (j3, Q) plane, pure hrt 
obviously corresponds to characteristics at constant temperature, and pure SCOZA 
to those at constant cutoff. We will therefore refer to these directions as HRT-like 
or SCOZA-Iike, respectively. 

For the pde (8), this directly implies the presence of a SCOZA-Iike characteristic 
for Q = Qoo and of an HRT-like one at f3 — 0; together they cover that part of 
the domain where the solution corresponds to the reference system. In addition, 
characteristics can have an HRT-like tangent only in points where ip(Q) = 0, and 
by considering all orders in the free parameter it is seen that all characteristics 
have HRT-like tangents whenever if)(Q) vanishes at Q < Qoo and j3 > 0. 

It is worth stressing that the significance of the characteristic curves extends 
not only to the numerics: In addition to the significance of shocks and rarefactions, 
neither of which are expected in the present setting, the pde is well-posed only 
when the characteristics establish a one-to-one mapping from the points at Q = 
to a subset of the boundary where the solution and its first derivatives are known. 
It is then possible to start out from the boundary and to integrate along those 
curves until the final solution is obtained at vanishing Q. 

But even for a stable and well-posed pde the characteristics are of vital im- 
portance for the numerics, independently of whether the discretization explicitly 
makes use of them: According to the Courant Friedrichs Levy criterion [17] any 
numerical scheme that does not adhere to the flow of information represented by 
the characteristics cannot be stable. 

4.2 Characteristics for the MSA-like case 

Even without prior knowledge of the generalized MSA solutions presented in sec- 
tion 3, some progress can be made towards the solution of eq. (8) by using the 
characteristics of the pde alone. In the MSA case the unknown parameter function 
X(x, y) of the preceding section corresponds to /i e ((3,Q). 
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4.2.1 Explicit results 

The situation is most transparent when considering only the limit of small z where 
eq. (8) becomes to order 0(z 2 ) 



dl_ 
3~ 



Q =F\Q)z, 



F(Q) =4nC [ 4>(k)k 2 dk. 

Jk>Q 



k>Q 

After elimination of z and a change of variables from Q to F we obtain a pde for 
I as a function of F and (3, viz., 

The characteristics of this linear first-order PDE are determined by 

dF _d/3 _ dJ 

which is trivially integrated to 

F(3 = const, , 21 , 
/ = const. 

The general solution of the pde is then any functional relation between the two 
constants of integration above, i. e., I can only depend on F(Q) j3, 

I = I(F(Q)j3). (22) 

This is also the case for the solutions of eq. (8) given in section 3 to first order in 
z. Eqs. (10) and (13) are both recovered to first order in z by imposing I = for 
P = 0. 

The situation is only slightly more complicated when terms of higher order 
in z are not neglected in eq. (8). We can obtain an equivalent PDE by cross- 
differentiating the equations (8) and eliminating d 2 I/dQd(3, as was done in sec- 
tion 4.1. From the relations (8) and the definition (11) we find 

-arc?. * q 7 L (HI =-JH) 

or 



i /an fdz\ i 

~z\dQ) z \Wq~P 



dJ_\ fdJ_\ fdz_\ 
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Rearranging this yields a quasilinear first-order pde for z(/3,Q), viz., 

1 fdJ_\ fdz_\ 1 fdJ\ ( dz_\ 1 (dJ_ 
~z \dQj, \W) Q ~ P \d~z) Q \dQj p ~P\dQ 

As in eq. (19), the equation for the characteristics becomes 

d(3 dQ dz 

1 fdJ_\ ~ ~ 1 7dT\ ~ 1 fd£\ ' 
z \dQj z p \J~z) Q \dQJ, 

This immediately gives zd(3 = /3dz so that 

fi e = - = const (24) 

along the characteristics. The relation between dQ and dz merely reproduces the 
total differential of J(Q,z), implying 

J = const. (25) 

The general solution of eq. (23) is then any functional relation between z//3 and 
J. It is easily seen that this is consistent with the expression for /x e given in 
section 3.2; constant /i e — z//3 and J = ji e K also imply constant K so that fi e 
must be a function of K only, just as indicated in eq. (14). 

4.2.2 Implications of the characteristics 

As pointed out in section 4.1, the pde has a well-defined and unique solution only 
when the characteristics establish a one-to-one mapping of the target system at 
various temperatures (Q — 0, (3 > 0) onto points where the solution is known 
from a suitable boundary condition. For the pde (8), the most natural condition 
to impose is provided by the reference system, the properties of which must be 
recovered both at f3 — and at Q = Qoo- On the other hand, among the general 
features of this type of pde in section 4.1 we found that these boundaries also 
coincide with the location of two characteristics. Consequently there can be no 
characteristic connecting a point on one of those boundaries with a target system 
at non-zero f3. At first sight this seems to render the reference system useless as 
a boundary condition. Indeed, for the MSA-like solution considered here, eq. (25) 
shows that the reference system only determines the solution where J = 0, i. e., 
on the boundary itself. 

However, this is not necessarily a severe problem as I {(3, Q) is expected to 
be continuous everywhere, including at the boundaries. The reference system 
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therefore still provides a valid description for systems removed from the boundaries 
by a very small amount, where the characteristic may have a direction different 
from that of the boundary and may actually lead to the target system at non-zero 
f3. In this case the reference system may still be used as a boundary condition; 
an example where even an infinitesimal separation is sufficient will be given in 
section 4.3. This situation echos the one found in continuum fluids in hrt: For 
those systems it has always been necessary to start numerical integration of the 
equations at some large but finite cutoff Qoo, typically on the order of 10 2 /cr, where 
a is a length characteristic of the repulsive hard core reference interaction [7, 18]. 
Neglecting the contribution of very high wavenumber fluctuations has never been 
a practical problem, nor is this expected to be the case here. 

In order to investigate whether a purely numerical solution is possible for the 
MSA-like equations considered in section 4.2.1, let us consider the results obtained 
by expanding the pde to first order in z, valid in the limit of high temperature. 
For the moment assuming ijj(0) ^ 0, there are potentials such that F(Q) vanishes 
only for Q = Q^. In this case / may be prescribed for all (3 at some fixed Q < Qoo, 
the characteristics stay at finite (3 at all lower Q, and the solution at Q = can 
be obtained at every temperature. 

On the other hand, F will be an oscillating function of Q for many short- 
ranged interactions, especially when the core condition is approximately taken 
into account by demanding ip(0) = 0, cf. section 2.2. At high temperature the 
field of characteristics then qualitatively looks as sketched in fig. 1. 

An important but unfortunate feature of that sketch are the SCOZA-Iike char- 
acteristics starting at f3 — whenever F(Q) = 0: These sectorals, as we shall call 
any characteristic the tangent of which at f3 = 0+ has a non-vanishing component 
in the (3 direction, divide the domain of the pde into a sequence of sectors. As 
there are no characteristics that cross from one sector into a neighbouring one, no 
information is transmitted across a sectoral, and the solutions in adjacent sectors 
are largely independent. The special significance of these characteristics comes 
from the fact that the solution for the target system (Q = 0) is determined by 
a boundary condition only when both lie in the same sector 2 . Furthermore, the 
SCOZA-like orientation of a sectoral means that the solution along it cannot depend 
on Fourier components of the interaction at other Q in the high temperature limit. 

For lower temperature we have to use the full solution. According to eq. (25), 
the sectorals are the curves where J = 0. As we further follow along them, their 
direction will no longer be SCOZA-like, and the solution along them will depend 
on a range of Fourier components. Still, one of the properties that can be inferred 

2 This is not the case for the analytical solution of section 3.2 as constancy of the A n is built 
right into eqs. (13) and (14). It is then possible to use a boundary condition to determine the 
A n in one sector, and to use their values in another one. 
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(3 = 0, reference system 




Figure 1: Sketch of the characteristics for the MSA-like solution to first order in 
z: Dashed lines show the characteristics whereas the dotted lines give the locus of 
■0 = where the characteristics have an HRT-like tangent. Sectorals are marked, 
as are the limits where the target and reference systems are recovered and the 
characteristics at f3 — and Q = Qoo- The inset gives the directions corresponding 
to pure hrt and SCOZA, respectively. There may or may not be a sectoral at Q — 0, 
depending on whether ijj(0) vanishes or not. As expected, there are no rarefactions, 
nor are there crossings of characteristics except at (3 — 0. 

from the temperature dependence of the J integral is that no sectoral can ever get 
to Q = at non-zero (3 since J > there. 

This feature of the sectorals also conforms to our expectation for a realistic 
closure relation. For suppose that some sectoral starts at (3 — and intermediate 
Q > and finally reaches Q = at (3 = (3 S > 0. Clearly it can only cover 
a finite Q interval so that the solution at Q = 0, (3 = (3 S cannot depend on 
Fourier components of the interaction outside this Q range. Absence of crossings 
of characteristics immediately implies that the same is also true for the properties 
of the target system at all (3 < (3 S . This situation is clearly unacceptable except 
for p a = 0. 

Another consequence of the temperature dependence of the J integral is the 
absence of sectorals at given Q for (3 above some Q-dependent limit; this limit 
is generally higher for larger Q. There are then only a few classes of admissible 
configurations of the characteristics in the ((3, Q) plane. Most likely, the sectorals 
are either driven towards infinite (3 at ever higher Q, or pairs of neighbouring 
sectorals join at their shared maximum (3 and so form loops; in the latter case, 
the maximum (3 occurs at a cutoff where ip(Q) = 0. Unfortunately, both of 
these possibilities are highly problematic: If sectorals stay separate but go towards 
infinite (3 and Q, there can be no finite temperature characteristic connecting a 
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point close to the reference boundary with one at Q = 0, (3 > 0. If, on the 
other hand, sectorals form loops, the same problem will arise at high cutoff. In 
addition, loops demarcate regions that are not entered by characteristics, implying 
rarefactions for (3 — > 0. 



4.3 Connecting reference and target systems 

The properties of the characteristics just inferred for the MSA case studied here are 
certainly unexpected and disappointing because they do not allow us to go from 
the reference system at Q — Qoo to the target system at vanishing Q. We now 
show that this is not so much a general defect of the combination of the energy 
and fluctuation routes but a consequence of the specific parameterization of the 
consistency problem. To see this, let us consider an arbitrary closure relation to 
the Ornstein-Zernike relation (3), giving the correlation functions c and g — h+ 1 
for any combination of Q, j3, and some unknown parameter A(/3, Q). For example, 
the MSA case considered so far corresponds to the direct correlation function 

c= --- + e(k-Q)v, A = ^. 
p n fi e 

Another simple possibility is 

£=--- + (A + Q(k - Q)) v. 
P P 

This is suggested by the usual hrt recipe for the continuum case when the core 
condition is not explicitly taken into account [4, 13, 19, 20]. 

In this more general situation with unspecified closure the correct pde is given 
by the earlier relations (4) and (5), except that the energy integral must be re- 
stricted to k > Q and that \x now depends not only on k, /3 and Q but also on the 
free parameter function X(/3,Q). As noted in section 2.2, the core condition must 
be fulfilled at every Q if the simple expression (5) is used, or else there will be an 
unphysical contribution to the internal energy from r = 0. 

In general, j2(k) is continuous at k = Q whereas the direct correlation function 
c(k) of the system corresponding to cutoff Q has a discontinuity of height v(Q) 
there. The relation between \x and c is obtained from eq. (6) by restricting the 
interaction to k < Q again, 

p + p 2 h(k) m 



i - fi{k) v{k) e(k - q) 



or 

m = --^rr:+v{k)Q{k-Q). 
P P\ k ) 
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Without further specifying the closure, a number of properties of the character- 
istics can easily be deduced: As the unknown parameter A enters eq. (8) only 
through c, we immediately see that the HRT-like characteristic at f3 — 0, the 
SCOZA-like characteristic at Q — Qoo, and the HRT-like direction of the tangents 
of all characteristics whenever if)(Q) = remain in this general setting. 

Both for discrete and continuous systems the question is then whether these 
boundary conditions at large but finite cutoff determine the solution at Q — 0. 
We have just seen that this is usually not the case for the pde (8). As pointed 
out before, the boundary conditions must be imposed in the same sector where 
the target system's properties are to be recovered. For this to be possible and 
the reference system to provide a valid initial condition for the integration of the 
characteristic equations, it is necessary but not sufficient that there be no sectorals 
above the initial value of Q, nor between the boundary condition and the target 
system. 

We then have to consider the existence and distribution of sectorals once again, 
which crucially depends on the precise way in which the free parameter enters 
the equations. As an example, let c ex (k;(3,Q) and h ex (k;(3,Q) be the direct and 
total correlation functions for any solution of the PDE such as, e. g., any of those 
presented in section 3. Inserting either of the two sample closures 

c(k; [3, Q, A) = c cx (k; P, Q) + A (fi(k) - fi(Q)) (26) 

and 

h(k; (3, Q, A) = h cx (k; (3,Q) + ^- (/ 2 (2fc -Q-Q^)- f 2 (Q + - 2k)) (27) 

(with largely arbitrary functions fi) into eqs. (4) and (5) yields a pde for \((3, Q). 
Obviously, the solution is A = for all (3 and Q in either case so that the original 
solution is reproduced, c = c cx and h = h ex . According to eq. (19), however, 
for eq. (26) the characteristics are HRT-like everywhere (d/3 = 0, so that there 
are no sectorals at all), whereas they are all SCOZA-like everywhere for eq. (27) 
(dQ = 0, so that there are infinitely many sectorals). Artificial and impractical 
as these examples are, they demonstrate that two different parameterizations of 
one and the same solution may have vastly different consequences, and eq. (26) in 
particular shows the potential for a closure free of sectorals. 

5 Necessary conditions for well-posedness 

We therefore conclude that the difficulties encountered in the case studied in sec- 
tion 4.2 are not inherent in the consistency problem but merely a consequence of 



17 



the particular MSA-like closure used. The results of our investigation thus leave 
open the possibility of a practical scheme, drawing our attention to the precise 
parameterization of the solution. 

While formulation of practical closure relations to be used together with eq. (8) 
lies outside the scope of the present contribution, we can give two general non- 
trivial conditions an ansatz for c must fulfill for the consistency problem to be 
well-posed and for a unique and differentiable solution to exist at all f3 > 0, 
Q > 0: (1) There must not be any sectorals in the sense of the definition given in 
section 4.3; and (2) any characteristic passing through a state with < (3 = 0(e) 
at some Q > must remain at a strictly positive inverse temperature (3 of order 
0(e) at lower cutoffs. While the former is quite obvious in the light of section 4.2.2, 
the second condition has not featured prominently so far: It merely expresses the 
absence of rarefactions and shocks at infinite temperature, and for the MSA case 
as treated here these already imply sectorals, cf. section 4.2.2. 

Both of these criteria share the advantage of involving only the limit (3 — > 0. 
They can therefore be checked by low order expansions in /3, which vastly simplifies 
analysis of the suitability of some specific closure relation in the context of eq. (8). 

Another important property is that both of them are quite general: Not only do 
they apply to discrete and continuous systems equally, they are also independent 
of the particular simplifications we chose to make in the present contribution. 
In particular, a closure may give a solution where the core condition is violated, 
leading to a spurious contribution to the energy integral of eq. (5). In this case 
more care should be exercised when evaluating the internal energy, leading to a 
modification of the pde [21]. The two conditions set out above, however, remain 
equally valid nevertheless. 

A careful study of the characteristics, especially of their compliance with the 
criteria just set out, is thus of vital importance for the consistency problem in- 
volving the energy and fluctuation routes only. The results of such an analysis are 
also highly relevant when the compressibility route enters the picture. For suppose 
that 1(13, Q, p) is the solution of the more elaborate pde so obtained. By evaluat- 
ing d 2 I /dp 2 and inserting this into the compressibility sum rule used in both pure 
HRT and pure SCOZA we arrive at a constraint involving both of the two unknown 
parameters in the parameterization of the solution at any ((3, Q, p). Restriction to 
fixed density in the spirit of the line method then again leads to a problem of the 
type considered here, and the properties of the characteristics of this restricted 
problem must be taken into account when solving, or discretizing, the PDE. In 
particular, the density dependence of the directions of the characteristics in the 
((3, Q) plane are then of prime importance for the numerical tractability of the pde 
by finite difference methods. At this point, however, it is not clear how to handle 
the situation where (3 does not increase monotonously along the characteristics, 
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nor do we know whether that case actually occurs with physically plausible closure 
relations and interactions. 

In summary, in the present contribution we have introduced and started to 
tackle the problem of unifying SCOZA and HRT, pointing out the importance of 
the condition of consistency between free and internal energies and for the time 
being setting aside consistency with the isothermal compressibility. A simplified 
lattice model system allowed us to derive a suitable pde and to study its known 
analytical solutions of MSA type. In general, however, we cannot assume an an- 
alytical solution to be available, and we therefore then turned to the problem of 
solving the pde without previous knowledge of the exact solution, especially with 
a view towards a numerical implementation: This highlighted the importance of 
the characteristic curves of the pde, and we have studied some of their properties 
and consequences. 
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